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HYPERCRITICAL ACCRETION ONTO A MAGNETIZED NEUTRON STAR 

SURFACE: A NUMERICAL APPROACH 

C. G. Bernal/ W. H. Lee/ and Dany Page^ 
RESUMEN 

La acrecion sobre una proto-estrella de neutrones en las horas que siguen 
al colapso del niicleo de una estrella masiva que le dio origen puede afectar 
sus propiedades observables. Este fenomeno se da en el regimen denominado 
hipcrcritico (Chevalier 1989), donde el enfriamiento por neutrinos es crucial para la 
evolucion termodinamica. En este trabajo presentamos un estudio en este contexto 
en una dimension con simetn'a esferica y llevamos a cabo simulaciones numericas en 
dos dimensiones dentro de una columna de acrecion sobre una estrella de neutrones. 
Consideramos procesos microfisicos detallados, enfriamiento por neutrinos y la pres- 
encia de campos magneticos en la aproximacion de magnetohidrodinamica ideal. 
Comparamos nuestros resultados numericos con las soluciones analiticas e inves- 
tigamos como las soluciones, tanto hidrodinamicas como magnetohidrodinamicas, 
difieren de estas. Iniciamos tambien una exploracion de como este proceso puede 
afectar la aparicion del remanente como un pulsar tipico en el radio. 

ABSTRACT 

The properties of a new-born neutron star, produced in a core-collapse su- 
pernova, can be strongly affected by the possible late fallback which occurs several 
hours after the explosion. This accretion occurs in the reg i me do minated by neu- 
trino cooling, explored initially in this context by [Chevalier Here we revisit 
this approach in a ID spherically symmetric model and carry out numerical sim- 
ulations in 2D in an accretion column onto a neutron star considering detailed 
microphysics, neutrino cooling and the presence of magnetic fields in ideal MHD. 
We compare our numerical results to the analytic solutions and explore how the 
purely hydrodynamical as well as the MHD solutions differ from them, and begin 
to explore how this may affect the appearance of the remnant as a typical radio 
pulsar. 

Key Words: Accretion — Hydrodynamics — Magnetic Fields — Supernovae: indi- 
vidual (SN1987A) — Stars: Neutron 



1. INTRODUCTION 



The neutrino signal (jHirata et al.l 119871: iBionta et al.l 119871) det ected from the supernova SN1987A clearly 
demonstrated the birth of a neutr on star (iBurrows k, Lattimeijri986[) . Identification of the pro genitor as the blue 
supergiant Sanduleak —69° 2 02a ( Gilmozzi et al. 1987t) and modeling of the early light curve ( Hillebrandt et al. 



1987 



Shigevama et al.lll987l ) proved that the supernova resulted from the core-collapse of a massive, ^ 20 Mg 



star. However, to date, there is no evidence for the pre sence of a pulsar, or even a quiet neutron star, in the 
remnant (see, e.g., discussion in lHaberl et al.ll2006l and Shternin fc Yakovlev 2008). Several soluti ons to this 
dilernma have been propose d as, e.g., the delayed collapse of th e neutron star into a black-hole (lEUis et al. 
19961: iBrown fc Bethel Il994l) or a delayed turn-on of the pulsar (iMichell [l993: iMuslimov fc Pagelll995l) . The 
latter case is just an extreme case of the more mundane possibility that the neutron star is weakly magnetized 
and/or slowly rotating, resulting in a spin-down energy that is low enough so as to be undetectable. 
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Fig. 1. Schematic evolution of a collapsing stellar core. A, the Chandrasekhar-mass Iron core collapses; B, upon 
formation of a proto-neutron star, the equation of state stiffens and produces a bounce, launching an outwardly prop- 
agating shock; C, a reverse shock is formed at the He-H interface in the stellar envelope, moving back towards the 
neutron star; D, an accretion shock is established, interior to which an atmosphere close to hydrostatic equilibrium in 
the neutrino-cooled regime deposits mass and energy onto the neutron star in the hours following the explosion. 



In recent years, several measurements of pulsar masses point toward large values approaching, or even 
exceeding, 2 Mq (see, e.g., [FroirG ,20081) which would strongly disfavor the black-hole explanation. On the 
other side, timing of radio quiet comp act stars in young supernova remnants, usually dubbed CCOs ("Central 
Comp act Objects" , Pavlov et al.ll2002 ) recently unveiled at least three case of weakly magnetized young neutron 
stars (jCotthelf fc HalpernI 120*081) : PSR J0821-4300 (in the SNR Puppis A) with a rotational period P = 112 
ms, an upper limit on its spi n-down power E < 2.3 x 10'^^ erg s^^, and a surface dipolar magnetic field 
strength B^p < 9.8 x 10" G (ICotthelf fc HalpernI l2009h : PSR 1E120 7 4-5209 (in the SNR PKS 1209-51/52) 
with P = 424 ms, E < 1.3 x lO^^ erg s'^, and B^ip < 3.3 x 10" G (jCotthelf fc HalpernI [20071) . and finally 



PSR J1852+0040 (i n the SNR Kes 79) with P = 424 ms, and measurements oi E = 3.0 x lO'^^ erg s'^, and 
Bdip = 3.1 X 10^° G (iHalpern fc Gottheljl2010l) . The last two of these hence have an energy output w ell below 



the 0.2-10.0 keV luminosity of the SN 1987A remnant, L < 5.7 x lO^* erg s'^ (jHaberl et al.H2006[) . If the 



neutron star produced by SN 1987A has similar characteristics it would presently be undetectable. 

In the present paper we consider the scenario in which the initial magnetic field of the new-born neutron 
star is str ongly modified by a phase of late, and intense, accretion, occurring a few hours after the initial 



explosion (iGeppert et al.lll999l). When a massive star explodes as a supernova, following the core-collapse 



scenario ( Wooslev fc Janka 2005 : Mezzacappa 2005t Janka et al. 2007), a large fraction of its mass expands 
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freely and interacts with the interstellar medium. However, the central compact remnant also interacts with 
the inner envelope through its gravitational field. In Type II SNe (see Figure [T]) , the initial core-collapse 
(panel A) produces a proto-neutron star when the equation of state stiffens close to nuclear density. A high- 
velocity (10^kms~^) expansive shock then starts moving outward (panel B). Flow lines bifurcate and some 
part of the matter falls back onto the central remnant. The rest is unbound and ends up being ejected in the 
explosion. This scenario produces a low density region in near free fall between the surface of the compact 
object and the extended atmosphere in near hydrostatic equilibrium that has formed. In case the progenitor 
star had a low density envelope surrounding the He core a reverse shock (panel C) decelerates the matter and 
causes a late fallback onto the compact object, depositing great amounts of matter onto t he surface of th e 
new-born neutr o n star in the hours followin g the explosion (panel D) . Followi ng the ideas of iBlondirJ (|l986f ). 
Chevalier ([i98i),|Houck k, Chevalieii (Il99ll) . andi rown &: Weingartneil (|l994l ) about the accretion of matter 



onto compact objects, it is possible to develop an analytical model of accretion following core-collapse, and 
particularly in the case of SN1987A. One of the salient features of this analysis is that the gas, being quite dense, 
is unable to cool by photon emission, and the mass accretion rates are highly super-Eddington in that sense. 
However, at sufficiently high temperatures, cooling through neutrinos sets in, mostly through pair annihilation 
and pair capture onto free nucleons, and given their much lower interaction cross section with matter, they 
are able to remove enough energy from the flow for accretion to take place. This regime is usually termed 
"hypercritical" accretion, and is common in the inner colla psing stellar cores and is likely to drive the central 
engines of Gamma Ray Bursts (|Lee fc Ramirez-Ruizll2007^ . With this model it is formally possible to obtain 
the radial position of the shock as a function of the mass accretion rate from fallback, assuming steady state 
in spherical symme try, a s well as the structure of t he envelope. 

Chevalier ( 1989f ) and Houck fc Chevalier ( 1991 ) computed such solutions in the context of SN1987A. Here 
we wish to explore the behavior of the flow under more general conditions, and present solutions for an 
accretion column in two dimensions, which we compare with the analytical scalings. Neutrino cooling is a 
crucial ingredient in the relevant density and temperature regimes, and we consider it along with a detailed 
equation of state. In addition and more importantly, we begin to explore the effects of the magnetic field on 
the accumulation of matter onto the neutron star surface. This is only possible through 2D simulations of the 
kind shown here, and we make a comparative analysis between the analytical and numerical approaches to 
consider the submergence of the magnetic field in the c rust of the neu t ron s tar and the piling up of matter 
on its surface. Previously, iMuslimov fc Page! (|l995[ ) and ICeppert et all (|l999l ) considered how such accretion 
might delay the switch-on of a pulsar following its formation in one-dim ensional cal c ulatio ns, computing the 
ohmic diffusion time of the magnetic field through the accreted matter. Frver et al. ( 1996[ ) studied the two- 
dimensional accretion dynamics onto new-born neutron stars in the neutrino cooled regime, finding that in some 
cases, neutrino-driven convection can significantly modify the simple one-dimensional steady state solution. 

Here we report on preliminary, two-dimensional numerical calculations which aim to determine if hypercrit- 
ical, neutrino-cooled accretion can submerge the magnetic field into the crust of the neutron star, and if it plays 
an important role in the dynamics in this regime. In § [2] we develop the analytical model of the hypercritical 
accretion process and calculate the structure of the envelope for a two dimensional accretion column. We build 
a numerical model, based on these analytical consideration, in § [3l In § [4l we show numerical results for various 
configurations including magnetic fields at several accretion rates and present a comparative analysis between 
the numerical and analytical solutions for the scenario of SN1987A. Finally, in §[5]wc present some preliminary 
conclusions. 



2. ANALYTICAL MODELS 

As a benchmark against which to compare our numer ical simula t ions, we summarize below the basic results 
of an analytical model, based on the one developed by Chevalier ( 19891 ). and adapt them to the case of an 
accretion column. The essential assumptions of the model are that the neutron star is at rest within the 
expanding medium at infinity and rotation is neglected. For this analytical approach, we also neglect the effect 
of a possible magnetic field. Matter is described by a polytropic equation of state, P — Kp'' with an index 
7 = 4/3, and is assumed to evolve adiabatically except at the shock interfaces and close to the neutron star 
surface where neutrino emission (through pair annihilation) assures that the accretion energy is properly 
removed from the system. 
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2.1. The Initial Late-Accretion Rate in SN1987A 



Spherically s ymmetric accre tion by a compact star in an initially static, infinite, background was described 
by iBondil (1952) (see also, e.g., IShapiro fc Teukolskv 1983 ). The mass accretion rate M is obtained from the 
density and sound velocity at infinity, poo and Coo respectively, as 



Mb = 47rA 



(1) 



where the numerical constant A = ~ 0.707 for the case of an ideal gas with adiabatic index 7 = 4/3. In 

our case the medium is not strictly initially at rest but has been set into expansion by the supernova shock 
wave. At early times the core is in homologous expansion, with a velocity v = r/t and density p such that 
pt'^ = Pata is constant, in terms of a reference density pa at time ta- As long as the time t is smaller than the 
Bondi accretion time scale T^ :^ GM/c^,, on e can still estimate M with Eq. ([T]) by allowing poo and 0^0 to be 



time dependent, and obtain (| Chevalier Il989l) 



^^B = 5.77^^(p,t3)V2t-3/2. 



According to IWooslevI (Il988l ) and IShigevama et al.l (|l988l) Pnf„ 10^ g 



(2) 



Now the densi t y and 



sound velocity at infinity have been estimated, for SN1987A, bv lWooslev ( 1988f ) and Bethe fc Pizzochero ( 1990[ ) 
as being of the order of 



2.5M 







1.78 X 10 



-13 



1.24 X 10'^ I — 



t 



gem 



9 —9 

cm s , 



(3) 
(4) 



where Vf ~ 600 kms"-'^ is the final expansion ve locity and Too ~ 70 keV(4 x 10^ cm/w/t)'^/'* for a radiation 
dominated shock as in a supernova like SN1987A. Bethe fc Pizzochero ( 1990f l calculated that the temperature 
for a shock radius of 4 x 10^ cm is Tsh — 70keV. The mass of the expanding CO core is ~ 4M0, but here we 
considered only 2.5 Mq because 1.5 M© were taken to make the compact object at the center of the supernova. 
We hence have 

15 15 

' gs-1 ~ 3.5 X 10-4 f:^) ' Mgyr-i. (5) 



M ~ 2.23 X 10^^ I — 



yr 



WooslevI (jl988l) calculated the time that the reverse shock takes to return to the surface of the neutron star for 
SN1987A as t ~ 7 X 10'^ s, giving us, for the accretion rate in SN1987A in the hypercritical regime 



M ~ 1.57 X lO^^gs"^ = 2500 Moyr^^ 



(6) 



This accretion rate exceeds by an order of magnitude the value calculated by iChevalierl (I1989D . M = 2.2 x 



^q28 g g-i _ 34Q yi~^ because of different assumed values at infinity. Now the Eddington mass accretion 
rate when considering photon radiation is MecU = 3.77 x lO^^gs"^, if one considers electron scattering in pure 
ionized Hydrogen as the source of opacit y, kes = .4cm ^ g~^. When M » Msdd the fiow is what we described 



above as Hypercritical Flow, studied by iBlondinI (|l986l ). For the case of SN1987A, we have 



M 



M. 



10^ - 10 



10 



(7) 



Edd 



for the two values given above, placing such flows clearly in the hypercritical, neutrino cooled regime. Henceforth 
we adopt as our fiducial accretion rate the value 



28 „, „-l 



Mo ~ 2.2 X lO^'^gs 



340 M© yr" 



(8) 
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2.2. The Envelope and the Shock Radius: Spherical Case 

When the reverse shock bounces against the surface of the neutron star, a third expansive shock is formed, 
which propagates through the infahing matter. Thus, eventuaUy an atmosphere in quasi-hydrostatic equihb- 
rium is formed around the compact object (see Panel D in Figure [I}, whose general structure can be calculated 
analytically under some simplif ying assumptions . 

Following the formulation of IChevaliediflOSgl) . the structure of the envelope is calculated and an expression 



for the pressure at the surface of the neutron star, P„s, in terms of M and the shock radius, r^/j, is derived. 
Cooling by neutrinos close to the neutron star surface, which depends on Pns , is introduced and this ultimately 
determines the shock radius solely as a function of the accretion rate. From the condition of hydrostatic 
equilibrium, dP/dr ~ —pGM/r'^, we obtain the integrated values for the pressure, density and velocity as 
a funtion of the distance from the neutron star surface. This is possible because neutrino cooling is only 
important near the surface of the neutron star and we can consider that the post-shock flow is adiabatic over 
the greater part of the volume. In addition, the flow is highly subsonic except close the shock. With this we 
obtain, P oc r^^/'^'i^^^ oc r^^,p oc r^'^/'^^^^^ cx r^^,v oc r'^^~'^^)/'^^~'^) (x r, where we have used 7 = 4/3. These 
results also are valid in the shock and the envelope structure in hydrostatic equilibrium is. 



r 



-4 



V = Vsh( — ), (11) 

where Psh, Psh, Vgh and Vsh are the values at the shock. These values can be obtained from the jump conditions 
when Psh >> Pq, where Pq and po the pre-shock pressure and density. Under these considerations we obtain 
Psh = {'^/8)poVo, Psh = 7po and Vsh = — (l/7)t'o- The pre-shock velocity is that of free-fall, vq = y^2GM/r^, 
and the density is po — M /ATTr^j^vo- From equation (8) we obtain the pressure at the surface, Pns = 1.36 x 
10^^^ Afr^^^, with M ~ 1.44M0 and r„s ^ 10^ c m. On the ot her hand, the energy loss by neutrinos (only pair 
production) by unit volume can be estimated as (lDicudll972l) . 



in = 1.83 X lO^-^'^P^-^^ergcm-^s-^ (12) 

In this case, we consider that e^ pairs also contribute to the pressure. Now, this cooling is operative only in a 
small volume close to the neutron star surface, ~ T^^ns since it is a sensitive function of temperature. So, from 
energy conservation, the shock radius is obtained as, 

10 

Tsh ^ 7.58 X 10^ — cm. (13) 

With this we have the structure of the envelope and the shock radius as a function of the accretion rate. 
For the case of SN1987A with our fiducial accretion rate M = Mq = 340 M© yr^^, the shock radius is 
rsh ^ 8.77 X 10^ cm. 



2.3. The Envelope and the Shock Radius: The Accretion Column 

If we consider a small rectangular accretion column of area A^oi onto a fraction of the neutron star surface, 
with area A = 47ri?^g, we can take it to be a plane-parallel surface (see Figure [2]). In this case, the spherical 
mass accretion rate must be scaled to its value in the column. Since in the spherical case the area depends on 
the distance to the neutron star, while for the case of an accretion column the area is constant, the structure 
of the envelope and the shock radius are modified. Note that this modification causes the mass accretion rate 
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Fig. 2. Schematic geometry of the accretion column onto the neutron star surface. 



per unit area to be independent of height in the domain. Also, since it is smaller than in the spherical case, 
the analytical estimate of the shock radius decreases significantly, and is now given by 

/ . X -0.16 

^/,, ^ 7.61 X loM ^'^[^ cm, (14) 

y Moyr V^coi ) 

where y measures the height above the neutron star surface. With these considerations and taking A^oi — 
(3 X 10^)^ cm^, the shock radius for the fiducial accretion rate, Eq. ([8|), is ysh — 6.92 x 10^ cm, and the 
structure of the envelope is, 

(15) 
(16) 
(17) 

The velocity profile is different for the accretion column as well because the column area is constant as a 
function of height above the neutron star. The conditions in the shock in the SN1987A scenario are thus 



= Psh 


- 




\ysh 






= Psh 






KVsh 






= Vsh 1 









vo = J ~ 7.53 X 10^cms-\ (18) 

M 

po ~ 2.31 X lO^gcm-3, (19) 

vq X A 

= 7po-1.62xl06gcm-3, (20) 

Vsh = -yi'o ^ -1.07 X 10^cms"\ (21) 

Psh = ^po«o- 1-14 X 10^5 dyncm-2. (22) 
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Now we can build a numerical model with more refined physics to perform 2D hydrodynamics (HD) and 
magnetohydrodynamics (MHD) simulations, which we can compare with the ID analytical results. 



3. NUMERICAL APPROACH 
For the work shown in this paper we used the numerical code AMR FLASH2.5 ( Frvxell et al. 2000[ ) to 



perform the 2D simulations. FLASH ( |http://fl ash.u chicago.edu/website/home/) is a modular, portable, highly 
scalable, adaptive-mesh simulation code for astrophysical hydrodynamics problems. It was originally developed 
at the DOE ASCI Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago for 
the purpose of simulating Type la supcrnovac, novae, and X-ray bursts. It has since evolved to handle more 
general astrophysical problems, including those involving coUisionless particle dynamics. FLASH is freely 
available from the ASCI Flash Center. This code is designed to allows users to configure initial and boundary 
conditions, change algorithms, and add new physics modules with minimal effort. It uses the PARAMESH 
library to manage a block-structured adaptative grid, placing resolution elements where they are needed most. 



3.1. The Numerical Method 



FLASH2.5 provides two main types of modules: Physics and Infrastructure Modules. In our model we used 
the hydro-mhd, eos-helmholtz, gravity and neutrino-cooling custom modules. 

The FLASH code solves the the equations of a magnetized fluid (ideal or non-ideal) , described by 



dp 

dt 



V-(pv) = 0, 



~dt 



V • (pvv - BB) + VP* = pg + V • T, 



dpE 

V- (v • r+crVT) + V • (B X (77V X B)) 
(9B 

V • (vB - Bv) 



V-[v(pP + P,)-B(vB)] = pvg+0(r,7/), 

= 0(T,r;), 
= -V X (77V X B) , 



where 



E 



dt 



1 



P2 

2p- 



= M 



(Vv) + (Vv)^ - ^ (Vv) 



(23) 

(24) 

(25) 
(26) 
(27) 

(28) 
(29) 
(30) 



Here P*, E and r are the total pressure, total specific energy and stress tensor, respectively, and the remaining 
symbols have their usual meaning. Units in these equations are such that no Air and /Uq factors appear. 

We have simplified the above set of equations by restricting ourselves to the ideal hydro and MHD cases. 
Setting the thermal conductivity, ct, and electrical resistivity, 77, to zero is justified by the fact that time 
scales for heat and magnetic field diffusion arc many orders of magnitude larger than our simulation times. 
The inviscid (p = 0) approximation, i.e., neglect of momentum diffusion, is acceptable because we have not 
considered rotation in our models. Note that when B = 0, the Euler equations are then obtained. 

A particular complication associated with solving the MHD equations numerically lies in the solenoidality 
of the magnetic field. The non-existence of magnetic monopoles, V • B = is difficult to satisfy in discrete 
computations. Being only an initial condition of the MHD equations, it enters the equations indirectly and is 
not, therefore, guaranteed to be generally satisfied unless special algorithmic provisions are made. FLASH2.5 
uses a simple yet very effective method to destroy the magnetic monopoles on the scale on which they are 
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generated. In this method, a diffusive operator proportional to VV • B is added to the induction equation, so 
that the equations become 

— + V • (vB - Bv) = -V X (r;V x B) - vV • B+Ty^ VV • B, (31) 

with the artificial diffusion coefHcient rja chosen to mach that of grid numerical diffusion. In the FLASH code, 
Va = (A/2)(1/Aa; + 1/Ay + 1/Az)^^, in 3D where A is the largest characteristic speed in the flow. Since 
the grid magnetic diffusion Reynolds number is always on the order of unity, this operator locally destroys 
magnetic monopoles at the rate which they are created. All our simulations are in cartesian coordinates: in the 
presence of a magnetic field polar/spherical coordinates are very troublesome and, presently, the MHD version 
of FLASH docs not support them. 



3.2. The Physics Ingredients 

In the analytical approach an ideal gas equation of state has been used. This allows much simplification in 
the structure of the envelope and in addition, the flow can be managed like an adiabatic fluid with 7 = 4/3. 
The gas is dominated by radiation, which is trapped within the flow. Also, the neutrino losses depend on a 
high power of the pressure, but are only important at the base of the envelope. Nevertheless, to account for 
the thermodynamics more accurately and for the consequent piling up of matter on the star, it is advisable and 
necessary to work with a more complete and realistic equation of state. The Hclmholtz EOS provided with 
the FLASH2.5 distribution contains more physics and is appropiate for addressing astrophysical phenomena in 
which electrons and positrons may be relativistic and/or degenerate and in which radiation may significantly 
contribute to the thermodynamic state. This EOS thus includes contributions from black-body radiation, 
completely ionized ideal nuclei, and free electrons and positrons. The pressure and internal energy are calculated 
as the sum over the components, 

Ptot — Prad Pion Pele ^~ Ppos Pcoul 
^tot — ^rad ~t~ ^ion ^~ ^ele ^pos ^~ ^coul- 

Here the subscripts "rad", "ion", "ele", "pos" and "coul" represent the contribution from radiation, nuclei, 
electrons, positrons, and Coulomb corrections, respectively. The radiation portion assumes a blackbody in local 
thermodynamic equilibrium, the ion portion (nuclei) is treated as an ideal gas with 7 = 5/3, and the electrons 
and positrons are treated as a non-interacting Fermi gas of arbitrary degeneracy and relativity. 

Under the physical conditions of interest for the set of simulations presented here, the gas is dense enough 
that the optical depth for photons is ^ 1, and they are fully trapped in the flow. Adding the corresponding 
term to the pressure as Prad = is thus entirely appropriate. We note that more recent versions of FLASH 

(upwards of 3.2) include modules for radiation transport, making them useful for a wider range of studies. 

The gravity module suplied with FLASH2.5 computes gravitational source terms for the code. These can 
take the form of the gravitational potential (/'(x) or the gravitational acceleration, 

g(x) = -V0(x). (34) 

The gravitational field can be externally imposed or self-consistently computed from the gas density via the 
Poisson equation, 

V^^(x) = 47rG'/9(x), (35) 

where G is Newton's gravitational constant. In the latter case, either periodic or isolated boundary conditions 
can be applied. In our case, we used an externally applied gravitational field (plane-parallel gravitational 
field), where the acceleration vector is parallel to one of the coordinate axes, and its magnitude drops with 
distance along that axis as the distance squared. Its magnitude and direction are independent of the other two 
coordinates. 

In the conditions present in both the high density part of the accretion flow and the underlying envelope 
neutrino emission occurs essentially through neutral currents processes. The five processes we included in the 
models are analogous to standard photon emission processes where the 7 emission is replaced by a — i7 pair. 
They are: 



(32) 
(33) 
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+ 17, 



the analogous of Compton scattering, 



PAIR ANNIHILATION: e" + e+ - 
PHOTONEUTRINOS: 7 + e± ^ f 
PLASMON DECAY: T ^ v + V, where P is a plasmon, 

BREMSSTRAHLUNG: + N ^ + N + v + V, where TV is a nucleus, and 
SYNCHROTRON: e± + B ^ + B + j/ + 77, where B represents the magnetic field. 
For t he first four processes we u sed the calculations of lltoh et alj (|l996r ) and for the synchrotron emission we 
followed iBezchastnov et al.l (|l997l ). Pair annihilation is the dominant process but synchrotron can make some 
significant contribution when the magnetic field becomes strongly compressed. As noted above, the density in 
the flow is typically high enough that photons are trapped, but not neutrinos. As a rough guide, the optical 
depth for neutrinos under coherent scattering ofi^ free nuclei is r^, ~ 1 when p ~ 10^^ g cm~^, which is several 
orders of magnitude higher than the maximum values studied here. Thus neutrino cooling can be implemented 
simply as a sink in the energy equation. 



3.3. The Initial and Boundary Conditions 

We simulated a small 2D accretion column in cartesian coordinates anchored onto the surface of the neutron 
star, and considered various accretion rates and magnetic field configurations. This set of simulations allows us 
to compare numerically obtained results in the pure hydrodynamical and MHD case with the proposed analytical 
approach, as well as to analyze the reaction of the magnetic field to the infalling gas. The computational domain 
covers < a; < 3 x 10^ cm, < y < 10^ cm. The dimensions for the column are: Acoi = Lx x L^ = (3 x 10^)^ cm^ 
(for the base) and heightcoi = Ly = 10^ cm for the height. This height is adequate because it is below the 
analytical shock radius value calculated for the accretion rate of SN1987A. The fluid is initially in free fall and 
we set a constant temperature in the gas. We considered horizontal {B^ ~ 10^^ G, By = 0), vertical {B^ = 0, 
By — 10^^ G, mimicking accretion onto the magnetic pole of the neutron star), diagonal {B^ = By = 10^^ G) and 
dipolar (Bx = 2fi/y^, By — 0, representing accretion onto the neutron star equator) cases, where /i = 5 x 10^^ 
is the dipolar moment, flxed so that B^ = 10^^ G at the neutron star surface. With these considerations, the 
initial conditions in the column for velocity, temperature and density are: 

P = — > (36) 

Vff X Acol 

T = 10^ K, (37) 

^ff = (38) 
B = 10^2 (39) 

where Mco\ — MAco\/A is the scaled accretion rate in the column making up the domain. 

For the vertical boundaries of the accretion column, x = and cc = 3 x 10^ cm, parallel to the y-axis, we 
implement standard periodic boundary conditions. Thus, any fluid element moving out of the computational 
domain on the right (left) boundary re-enters the domain on the left (right) edge with the same thermodynamical 
properties and velocity. 

For the top and bottom of the computational domain, parallel to the x-axis, we implemented custom 
boundary conditions. The gravity vector is along the y-axis, and we want the lower boundary at y = to 
support the fluid above against infall, mimicking the hard surface of the neutron star, and in addition to have 
the magnetic fleld anchored to it. In order to establish this boundary, we use "guard", or "ghost" numerical 
cells. These are cells outside the formal computational domain (e.g., at y < or y > 10^) for which we can flx 
the hydrodynamical and thermodynamical properties and that are not evolved along with the rest of the flow. 
They are useful precisely to guarantee boundary conditions of interest, depending on the setup of the problem. 
A layer of at least 2 such cells along the top and bottom of the domain can thus be used to compute proper 
gradients at the edge of the flow (e.g., a pressure or temperature gradient). In this case, along the bottom 
edge of the column, we flx the velocities to be null in all guard cells, {vx = Vy = Vz = 0)-, keeping them at rest, 
and copy the density and the pressure of the flrst zone of the numerical domain to mimick the the neutron 
star surface: (p = /c(l), P = P{^) + pv'^ + pgh), where the label 1 refers to the flrst cell in the computational 
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domain. The magnetic field is put in this boundary in such form that it is continuous from the guard cell to 
the physical domain, i.e, we anchor the magnetic field onto the neutron star surface and in the rest of the guard 
cells it is null. The other thermodynamics variables are calculated from the equation of state. At the top of 
the column, y ~ 10^ cm, we set the velocity to be that of free fall, Vy = — ^ 2GM j y in all the guard cells, and 
set the density to fix a constant inflow mass accretion rate, p = Mcovj {\vy\ x Acoi). As in the computational 
domain initially, the temperature in the guard cells is set to T = 10^ K (at all times). The remaining variables 
are calculated from the equation of state. 

4. RESULTS AND DISCUSSION 

We now present results obtained from the 2D hydrodynamical simulations (HYDRO) as well as for the 
MHD case for an accretion column in cartesian coordinates, and compare these to the analytical scalings. We 
varied the accretion rate and magnetic field configuration (for the MHD case). The chosen rates were one, two 
and three orders of magnitude above our fiducial rate Mq — 340 Mq yr~^. 

4.1. Comparison of the HYDRO and MHD solvers 

For the assumed physical parameters of SN1987A, in 200 ms the system reaches a quasi-stationary state, 
whereas for higher rates of accretion, this drops substantially: 60 ms at IOMqi 20 ms at IOOMq and 5 ms at 
lOOOMg- We set a level of refinement of 4, with 2 blocks along the x-axis and 18 along y-axis. This implies an 
effective resolution of 128 x 1152 zones in the computational domain. 

In Fig. Owe show the density contrast for the HYDRO and MHD cases with null magnetic field (MHD_0), 
for M = lOOAfo- We choose this accretion rate as being representative since its associated shock radius is 
much smaller than for M = Mq, and is therefore easier to visualize. In addition, it is possible to both do a 
comparative analysis of solvers (HYDRO and MHD) and of their response to the imposed initial conditions. 
This comparison allows us to determine whether the equations are being solved in both modules to a comparable 
accuracy. In principle, the MHD module with null magnetic field should reproduce exactly the results obtained 
with module HYDRO. The constrasts of pressure, specific total energy, velocity and neutrino cooling per unit 
volume for all the cases (at t = 10 ms), are shown in Fig. The radial profiles of density, pressure and velocity 
for the SN1987A accretion rate are given in Fig. [5] 

We note that although the system reaches a quasi-stationary state in t = 20 ms, there is remnant noise 
in the radial profile of the velocity due to the interaction of the matter with the lower boundary condition 
and to the fact that horizontal motions are allowed because of the periodic boundary condition. On the other 
hand, only the bottom section, 4 x 10^ cm, of the entire accretion column, with height lO'^ cm is shown, where 
the most interesting processes occur. We note that the profiles, while not identical in all respects, are indeed 
very similar, showing that the HYDRO and MHD solvers are giving essentially the same final state, both in 
space and time evolution. There is some convection in the early stages of the evolution, and Rayleigh- Taylor 
instabilities are present, but are quickly damped as the system approaches the stationary solution. Deviations 
from this are most evident when one examines variations in the velocity field. 

The location of the shock is reproduced quite well, to within 5% when compared to the analytical calculation. 
Moreover, both solvers place it essentially at the same height, indicating that the quantitative aspects are not 
affected from one to the other. Since the code is able to model the bottom of the column self-consistently within 
the imposed boundary condition, the numerical solution deviates from the self-similar scaling once neutrino 
cooling becomes important, and matter starts piling up near the surface. 

Hereafter, unless otherwise noted we refer to calculations with M = Mq. The adiabatic and radiative 
gradients can be calculated from the simulations, when the system is relaxed. We find Vad = 1 — l/7c — 0.26, 
Vrad = (dlnT/dlnP) ~ 0.24. In this case, the value of the adiabatic index 7c has been taken directly from 
the simulation (7c = 1.35), and the radiative gradient was calculated by building a plot of temperature vs. 
pressure. These gradients have almost constant values within the envelope, except in the region close to the 
neutron star surface. Since Vad > Vrad, the system is manifestly stable to convection. Nevertheless, being so 
close numerically is probably indicative of marginal stability. Within the envelope the flow is fully subsonic, as 
expected after passing through the accretion shock front: the sound speed is Cg = \/lcP/p — 6.9 x 10^ cms^^, 
and V ~ 1.24 x lO^cms"^, giving a Mach number m = v/cg ~ 10^^. Therefore, besides confirming that the 
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HYDRO and MHD solvers give accurate and consistent results, we are able to study the global structure of 
the accretion column in detail and compare it with the analytical approach, particularly in the region where 
the approximations in the latter break down. 

It is worthy to note the thermodynamical conditions the fluid is in as it accretes towards the proto-neutron 
star. The Fermi temperature can be computed from the Fermi energy Ep 



Tp^ — = — - — — ~ 6.48 X 10^°K, 

KB KB 
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Fig. 4. Color maps of pressure, total specific energy, magnitude of velocity and neutrino emissivity, for cases HYDRO 
(top) and MHD_0 (bottom) at t = 10 ms. The accretion rate is M = lOOMo- 



at the base of the flow, where pp = (37r^ne)^/"^?i is the Fermi momentum. The temperature obtained from 
the simulation, close to the bottom of the aceretion column in quasi-stationary state is T ~ 4.54 x 10^" K, 
so T/Tp ~ 0.7. It is thus clear that assuming that the e"^ pairs are entirely degenerate is not a proper 
approximation, and a full expression such as the one in the Helmholtz equation of state is required if one 
wishes to compute the evolution of the flow accurately. It is also clear that neutrino cooling effectively turns 
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Fig. 5. Radial profiles of density (top), pressure (middle) and magnitude of velocity (bottom) at t = 200 ms. The 
accretion rate is M = Mo. The HYDRO (red), MHD_0 (green) and analytical solution (blue) are shown together. The 
location of the shock is reproduced with very good agreement in the two numerical cases with respect to the analytical 
solution. In the shocked region, the velocity shows significantly higher behavior due in part to lateral motions of the 
gas. At small radii, the numerical solutions deviate from the analytical curve, since the assumption of self-similarity is 
no longer valid as material piles up near the star and cooling becomes relevant in the energy balance. 



on at a scale height Ly ^ 2 x 10^ cm. For the simulation with AI = A/q, the integrated neutrino lumi nosity, 
shown in Fig. [6l is L^, ~ 2.51 x lO'^^ergs"^, close to the value estimated with the cooling function of iDicus 
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Fig. 6. Neutrino luminosity integrated over the computational domain for the fiducial accretion rate, M = Mo, up to 
t — 200 ms. Note how after an initial transient, the power output is leveling off as the system reaches a quasi-stationary 
state. 



(|l972l ) scaled to the column: = £„ x F ~ 1.83 x 10''6ergs-\ with V ~ {2 x 10^) x (3 x 105)2cm3. 

Once the system reaches the quasi-stationary state, radial profiles can be compared for different accretion 
rates. Four different rates for each initial condition were computed. In all of these, the piling up of material close 
to the neutron star surface is seen. The velocity profiles remain noisy and turbulent in the shocked region, but on 
average the analytical profile is globally recovered. In Fig. [7]these are plotted, along with density and pressure, 
for case MHD_0. Note also that at greater accretion rates the shock is located at lower height, as expected. For 
M/Mo = 1,10,100,1000, the position of the shock in the simulation is at Rsh/W^cm = 7.06,4.96,3.74,2.68, 
in excellent agreement with the analytical values given by i?s/j/10^cm = 6.92, 4.80, 3.33, 2.31, respectively (see 
Fig. El). 



4.2. Magnetic field submergence 

We now consider the case with non-zero magnetic field strength. Fig. |9] shows the radial profiles of density, 
pressure and magnitude of the velocity for M = Mq with several field configurations: null (MHD_0), constant 
horizontal (MHD_H), constant vertical (MHD_V), constant diagonal (MHD_D) and dipolar (MHD.DIP), for 
comparative effects. The initial intensity of the magnetic field in all cases is 10^^ G, except in the dipole 
configuration, where it is 10^^ G at the neutron star surface. We also overplot the hydrodynamical solution for 
comparison. Note that the profiles are practically the same at this accretion rate indicating that the magnetic 
field is not playing an important role as far as the dynamics are concerned. 

In all simulated cases, regardless of the magnetic field configuration, when the system has relaxed and 
reached the quasi-stationary state, the field is completely submerged in the neutron star crust. Its intensity 
rises accordingly, by up to two orders of magnitude for the highest accretion rates. Fig. llOl shows the distribution 
of magnetic field strength after the system has relaxed, when M = 1000 Mq, for our four initial magnetic field 
configurations. It is only within the first km in the column, where the matter piles up, that the magnetic field 
is at or above the initial value in the calculation, and the compression is quite clear. 

The initial dynamics in the MHD case are somewhat more violent than in the pure hydrodynamical case. 
The infalling gas quickly drags the initial field towards the neutron star surface since the ram pressure. Pram = 
pv'^/2 is substantially greater than the magnetic pressure Pmag = B'^/Stt, even for the smallest accretion rate, 
M = Mq. The increased magnetic pressure as compression takes place is insufficient to overcome this flow, 
and large field strengths close to the surface result. The effect on the large scale dynamics is thus of a more 
transitory nature, and sensitive to the initial conditions, than a permanent feature. As a second point, we note 
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Fig. 7. Radial profiles of density (top), pressure (middle) and velocity (bottom) for the MHD_0 case after a stationary 
state has been reached. The accretion rates are M / Mq = 1, 10, 100, 1000 (cyan, blue, green and red, respectively). 



that the magnetic field, advected along with the flow, fluctuates in strength strongly in the shocked region as 
it piles up against the lower boundary, where neutrino cooling is eflicient. The additional piling up of material 
makes it even harder for the fleld to rise to significant levels as the evolution proceeds further. Nevertheless, as 
the system evolves the turbulent structures that form initially begin to smooth themselves until they disappear 
completely in the hydrodynamical case, but some small scale structure remains when magnetic fields are present. 
Once the accretion rate drops significantly, it is in principle possible that the field will rise buoyantly through 
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Fig. 8. Accretion rate as a function of shock radius for the hydrodynamical simulations with M /Mo ~ 1000, 100, 10, 1 
after a stationary state has been reached. The analytical solution given by equation [T?] is shown for comparison. 



the envelope, playing some dynamical role as the accretion time becomes long and the balance between ram 
and magnetic pressure is reversed. This will occur on a much longer time scale than simulated here, and its 
modeling requires a different set of assumptions in terms of the present set of calculations. 



5. CONCLUSIONS 



We have presented the results of two-dimensional simulations of accretion in the hypercritical, neutrino- 
cooled regime onto the surface of a neutron star, using the FLASH code. The flow in accretion columns for a 
variety of initial accretion rates was simulated until a steady state was reached. We find that at this stage, the 
location of the accretion shock, where the flow transitions from free fall to subso nic settling onto the neutron 
star surface, is well reproduced when compared with the analytical estimates of IChevalier (Il989l) . However, 
close to the surface, matter piles up, the solution is no longer adiabatic, and the self-similar character of the 
flow breaks down as expected. 

We performed a detailed comparison of the hydrodynamical and ideal MHD routines in FLASH, and found 
excellent agreement between the two when the initial field is null. For various finite field configurations (initially 
horizontal, vertical, diagonal and dipolar), we find that performing the calculations in two dimensions does 
not allow for any additional buoyancy effects of the field to be manifested: for all accretion rates simulated, 
the initial field is entirely advectcd by the fiow and submerged close to the neutron star surface. Its intensity 
rises accordingly, by up to two orders of magnitude in some cases. In principle, thus, it is possible for such 
an accretion episode following core collapse and the formation o f a proto-neutron stars t o effectively bury the 
initial field and delay the appearance of a classical radio pulsar (jMuslimov &: Pagelll995l ). The simulated time 
scales at present do not allow us to place hard constraints on the re-diffusion of the field at late times, and a 
more quantitative estimation of this is left for future work. 
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